#Follow the Fox: Cross-sectional coding, data analysis, and robustness checks

#######################
#Packages

library(tidyverse)
library(haven)
library(dotwhisker)
library(jtools)
library(sjPlot)
library(DataExplorer)
library(psych)
library(texreg)
library(ggeffects)
library(ggplot2)
library(forcats)
library(ltm)
library(lavaan)
library(hbal)
library(sensemakr)

# load data
setwd("#Your_Working_Directory")

load("FollowFox_R_Cross_Sectional_Wave3.Rdata")

load("FollowFox_R_Cross_Sectional_Wave4.Rdata")

#recode Wave 3 data
wave3 <- wave3%>%
  mutate(
    ABC = case_when(
      Q11_m_W3_1 == 1 ~ 1,
      Q11_m_W3_1 == 2 ~ 0,
      Q11_m_W3_1 == 9 ~ NA_real_),
    CBS = case_when(
      Q11_m_W3_2 == 1 ~ 1,
      Q11_m_W3_2 == 2 ~ 0,
      Q11_m_W3_2 == 9 ~ NA_real_),
    NBC = case_when(
      Q11_m_W3_3 == 1 ~ 1,
      Q11_m_W3_3 == 2 ~ 0,
      Q11_m_W3_3 == 9 ~ NA_real_),
    CNN = case_when(
      Q11_m_W3_4 == 1 ~ 1,
      Q11_m_W3_4 == 2 ~ 0,
      Q11_m_W3_4 == 9 ~ NA_real_),
    Fox_News = case_when(
      Q11_m_W3_5 == 1 ~ 1,
      Q11_m_W3_5 == 2 ~ 0,
      Q11_m_W3_5 == 9 ~ NA_real_),
    MSNBC = case_when(
      Q11_m_W3_6 == 1 ~ 1,
      Q11_m_W3_6 == 2 ~ 0,
      Q11_m_W3_6 == 9 ~ NA_real_),
    PBS = case_when(
      Q11_m_W3_7 == 1 ~ 1,
      Q11_m_W3_7 == 2 ~ 0,
      Q11_m_W3_7 == 9 ~ NA_real_),
    Local_News = case_when(
      Q11_m_W3_8 == 1 ~ 1,
      Q11_m_W3_8 == 2 ~ 0,
      Q11_m_W3_8 == 9 ~ NA_real_),
    Spanish_News = case_when(
      Q11_m_W3_9 == 1 ~ 1,
      Q11_m_W3_9 == 2 ~ 0,
      Q11_m_W3_9 == 9 ~ NA_real_),
    Newsmax = case_when(
      Q11_m_W3_10 == 1 ~ 1,
      Q11_m_W3_10 == 2 ~ 0,
      Q11_m_W3_10 == NA ~ 0),
    No_News = case_when(
      Q11_m_W3_11 == 1 ~ 1,
      Q11_m_W3_11 == 2 ~ 0,
      Q11_m_W3_11 == 9 ~ NA_real_),
    GRT1 = case_when(# 5: Strongly disagree, 1: Strongly agree
      Q41_1_W3 == 5 ~ 0,
      Q41_1_W3 == 4 ~ .25,
      Q41_1_W3 == 3 ~ .5,
      Q41_1_W3 == 2 ~ .75,
      Q41_1_W3 == 1 ~ 1,
      TRUE ~ NA_real_),
    GRT2 = case_when(
      Q41_7_W3 == 5 ~ 0,
      Q41_7_W3 == 4 ~ .25,
      Q41_7_W3 == 3 ~ .5,
      Q41_7_W3 == 2 ~ .75,
      Q41_7_W3 == 1 ~ 1,
      TRUE ~ NA_real_),
    GRT3 = case_when(
      Q41_10_W3 == 5 ~ 0,
      Q41_10_W3 == 4 ~ .25,
      Q41_10_W3 == 3 ~ .5,
      Q41_10_W3 == 2 ~ .75,
      Q41_10_W3 == 1 ~ 1,
      TRUE ~ NA_real_),
    GRT_Beliefs = (GRT1+GRT2+GRT3)/3,
    Male = case_when(# 1: Male, 2: Female
      gender_W3 == 1 ~ 1,
      gender_W3 == 2 ~ 0),
    White = if_else(raceOS_W3==1,1,0),
    Education = case_when(# 1: least educ, 6: most educ
      educ_W3 == 1 ~ 0,
      educ_W3 == 2 ~ .2,
      educ_W3 == 3 ~ .4,
      educ_W3 == 4 ~ .6,
      educ_W3 == 5 ~ .8,
      educ_W3 == 6 ~ 1),
    Income = case_when(# 1: least, 16: most
      faminc_new_W3 == 1 ~.0625,
      faminc_new_W3 == 2 ~.125, 
      faminc_new_W3 == 3 ~.1875,
      faminc_new_W3 == 4 ~.25, 
      faminc_new_W3 == 5 ~.3125, 
      faminc_new_W3 == 6 ~.375,
      faminc_new_W3 == 7 ~.4375,
      faminc_new_W3 == 8 ~ .50,
      faminc_new_W3 == 9 ~.5625, 
      faminc_new_W3 == 10 ~.625,
      faminc_new_W3 == 11 ~.6875,
      faminc_new_W3 == 12 ~.75,
      faminc_new_W3 == 13 ~.8125,
      faminc_new_W3 == 14 ~.875,
      faminc_new_W3 == 15 ~.9375,
      faminc_new_W3 == 16 ~1,
      TRUE ~ NA),
    Ideology = case_when(# 1: liberal, 5: conservative
      ideo5_W3 == 1 ~ 0,
      ideo5_W3 == 2 ~ .25,
      ideo5_W3 == 3 ~ .5,
      ideo5_W3 == 4 ~ .75,
      ideo5_W3 == 5 ~ 1,
      TRUE ~ NA_real_),
    Party_ID = case_when(# 1: Democrat, 7: Republican
      pid7_W3 == 1 ~ 0,
      pid7_W3 == 2 ~ .1666667,
      pid7_W3 == 3 ~ .3333333,
      pid7_W3 == 4 ~ .5,
      pid7_W3 == 5 ~ .6666667,
      pid7_W3 == 6 ~ .8333333,
      pid7_W3 == 7 ~ 1,
      TRUE ~ NA_real_),
    Age = (2024-birthyr_W3)/100,
    Trump_Therm = Q62_2_W3/100,
    Republican_Therm = Q34_6_W3/100,
    FIRE1 = case_when(Q53_4_W3==1~0,
                      Q53_4_W3==2~.25,
                      Q53_4_W3==3~.5,
                      Q53_4_W3==4~.75,
                      Q53_4_W3==5~1),
    FIRE2 = case_when(Q53_5_W3==1~1,
                      Q53_5_W3==2~.75,
                      Q53_5_W3==3~.5,
                      Q53_5_W3==4~.25,
                      Q53_5_W3==5~0),
    FIRE3 = case_when(Q53_6_W3==1~0,
                      Q53_6_W3==2~.25,
                      Q53_6_W3==3~.5,
                      Q53_6_W3==4~.75,
                      Q53_6_W3==5~1),
    FIRE_Scale = (FIRE1+FIRE2+FIRE3)/3)

wave3 <- wave3 %>%
  mutate(ABC = if_else(is.na(ABC), 0, ABC),
         NBC = if_else(is.na(NBC), 0, NBC),
         CBS = if_else(is.na(CBS), 0, CBS),
         Fox_News = if_else(is.na(Fox_News), 0, Fox_News),
         CNN = if_else(is.na(CNN), 0, CNN),
         MSNBC = if_else(is.na(MSNBC), 0, MSNBC), 
         PBS = if_else(is.na(PBS), 0, PBS),
         Local_News = if_else(is.na(Local_News), 0, Local_News),
         Spanish_News = if_else(is.na(Spanish_News), 0, Spanish_News),
         Newsmax = if_else(is.na(Newsmax), 0, Newsmax), 
         No_News = if_else(is.na(No_News), 0, No_News))

#add recoded variables for continuous Fox Measure
wave3 <- wave3 %>%
  mutate(Fox1 = case_when(Q22_m_W3_1 ==1 ~ 1,
                          Q22_m_W3_1 ==2 ~ 0,
                          TRUE ~ 0),
         Fox2 = case_when(Q22_m_W3_2 ==1 ~ 1,
                          Q22_m_W3_2 ==2 ~ 0,
                          TRUE ~ 0),
         Fox3 = case_when(Q22_m_W3_3 ==1 ~ 1,
                          Q22_m_W3_3 ==2 ~ 0,
                          TRUE ~ 0),
         Fox4 = case_when(Q22_m_W3_4 ==1 ~ 1,
                          Q22_m_W3_4 ==2 ~ 0,
                          TRUE ~ 0),
         Fox5 = case_when(Q22_m_W3_5 ==1 ~ 1,
                          Q22_m_W3_5 ==2 ~ 0,
                          TRUE ~ 0),
         Fox6 = case_when(Q22_m_W3_6 ==1 ~ 1,
                          Q22_m_W3_6 ==2 ~ 0,
                          TRUE ~ 0),
         Fox7 = case_when(Q22_m_W3_7 ==1 ~ 1,
                          Q22_m_W3_7 ==2 ~ 0,
                          TRUE ~ 0),
         Fox8 = case_when(Q22_m_W3_8 ==1 ~ 1,
                          Q22_m_W3_8 ==2 ~ 0,
                          TRUE ~ 0),
         Fox9 = case_when(Q22_m_W3_9 ==1 ~ 1,
                          Q22_m_W3_9 ==2 ~ 0,
                          TRUE ~ 0),
         Fox10 = case_when(Q22_m_W3_10 ==1 ~ 1,
                           Q22_m_W3_10 ==2 ~ 0,
                           TRUE ~ 0),
         Fox_Total = (Fox1+Fox2+Fox3+Fox4+Fox5+Fox6+Fox7+Fox8+Fox9+Fox10)/10)

#add data on political violence
wave3 <- wave3 %>%
  mutate(TrumpVoter = case_when(Q67_W3==1~1,
                                Q67_W3==2~0,
                                TRUE~NA_real_),
         Violence1 = case_when(Q64b_1_W3==1~1,
                               Q64b_1_W3==2~0.75,
                               Q64b_1_W3==3~0.5,
                               Q64b_1_W3==4~0.25,
                               Q64b_1_W3==5~0,
                               Q64b_1_W3==6~NA_real_),
         Violence2 = case_when(Q64b_2_W3==1~1,
                               Q64b_2_W3==2~0.75,
                               Q64b_2_W3==3~0.5,
                               Q64b_2_W3==4~0.25,
                               Q64b_2_W3==5~0,
                               Q64b_2_W3==6~NA_real_),
         Violence3 = case_when(Q64b_3_W3==1~1,
                               Q64b_3_W3==2~0.75,
                               Q64b_3_W3==3~0.5,
                               Q64b_3_W3==4~0.25,
                               Q64b_3_W3==5~0,
                               Q64b_3_W3==6~NA_real_))

#add data on group attitudes
wave3 <- wave3 %>%
  mutate(Latino_Therm = Q34_1_W3/100,
         White_Therm = Q34_3_W3/100,
         White_Preference = White_Therm - Latino_Therm)

#Recode Wave 4 data

wave4 <- wave4%>%
  mutate(
    ABC = case_when(
      Q11_m_W4_1 == 1 ~ 1,
      Q11_m_W4_1 == 2 ~ 0,
      Q11_m_W4_1 == 9 ~ NA_real_),
    CBS = case_when(
      Q11_m_W4_2 == 1 ~ 1,
      Q11_m_W4_2 == 2 ~ 0,
      Q11_m_W4_2 == 9 ~ NA_real_),
    NBC = case_when(
      Q11_m_W4_3 == 1 ~ 1,
      Q11_m_W4_3 == 2 ~ 0,
      Q11_m_W4_3 == 9 ~ NA_real_),
    CNN = case_when(
      Q11_m_W4_4 == 1 ~ 1,
      Q11_m_W4_4 == 2 ~ 0,
      Q11_m_W4_4 == 9 ~ NA_real_),
    Fox_News = case_when(
      Q11_m_W4_5 == 1 ~ 1,
      Q11_m_W4_5 == 2 ~ 0,
      Q11_m_W4_5 == 9 ~ NA_real_),
    MSNBC = case_when(
      Q11_m_W4_6 == 1 ~ 1,
      Q11_m_W4_6 == 2 ~ 0,
      Q11_m_W4_6 == 9 ~ NA_real_),
    PBS = case_when(
      Q11_m_W4_7 == 1 ~ 1,
      Q11_m_W4_7 == 2 ~ 0,
      Q11_m_W4_7 == 9 ~ NA_real_),
    Local_News = case_when(
      Q11_m_W4_8 == 1 ~ 1,
      Q11_m_W4_8 == 2 ~ 0,
      Q11_m_W4_8 == 9 ~ NA_real_),
    Spanish_News = case_when(
      Q11_m_W4_9 == 1 ~ 1,
      Q11_m_W4_9 == 2 ~ 0,
      Q11_m_W4_9 == 9 ~ NA_real_),
    Newsmax = case_when(
      Q11_m_W4_10 == 1 ~ 1,
      Q11_m_W4_10 == 2 ~ 0,
      Q11_m_W4_10 == NA ~ 0),
    No_News = case_when(
      Q11_m_W4_11 == 1 ~ 1,
      Q11_m_W4_11 == 2 ~ 0,
      Q11_m_W4_11 == 9 ~ NA_real_),
    GRT1 = case_when(# 5: Strongly disagree, 1: Strongly agree
      Q41_1_W4 == 5 ~ 0,
      Q41_1_W4 == 4 ~ .25,
      Q41_1_W4 == 3 ~ .5,
      Q41_1_W4 == 2 ~ .75,
      Q41_1_W4 == 1 ~ 1,
      TRUE ~ NA_real_),
    GRT2 = case_when(
      Q41_7_W4 == 5 ~ 0,
      Q41_7_W4 == 4 ~ .25,
      Q41_7_W4 == 3 ~ .5,
      Q41_7_W4 == 2 ~ .75,
      Q41_7_W4 == 1 ~ 1,
      TRUE ~ NA_real_),
    GRT3 = case_when(
      Q41_10_W4 == 5 ~ 0,
      Q41_10_W4 == 4 ~ .25,
      Q41_10_W4 == 3 ~ .5,
      Q41_10_W4 == 2 ~ .75,
      Q41_10_W4 == 1 ~ 1,
      TRUE ~ NA_real_),
    GRT_Beliefs = (GRT1+GRT2+GRT3)/3,
    Male = case_when(# 1: Male, 2: Female
      gender_W4 == 1 ~ 1,
      gender_W4 == 2 ~ 0),
    White = if_else(raceOS_W4==1,1,0),
    Education = case_when(# 1: least educ, 6: most educ
      educ_W4 == 1 ~ 0,
      educ_W4 == 2 ~ .2,
      educ_W4 == 3 ~ .4,
      educ_W4 == 4 ~ .6,
      educ_W4 == 5 ~ .8,
      educ_W4 == 6 ~ 1),
    Income = case_when(# 1: least, 16: most
      faminc_new_W4 == 1 ~.0625,
      faminc_new_W4 == 2 ~.125, 
      faminc_new_W4 == 3 ~.1875,
      faminc_new_W4 == 4 ~.25, 
      faminc_new_W4 == 5 ~.3125, 
      faminc_new_W4 == 6 ~.375,
      faminc_new_W4 == 7 ~.4375,
      faminc_new_W4 == 8 ~ .50,
      faminc_new_W4 == 9 ~.5625, 
      faminc_new_W4 == 10 ~.625,
      faminc_new_W4 == 11 ~.6875,
      faminc_new_W4 == 12 ~.75,
      faminc_new_W4 == 13 ~.8125,
      faminc_new_W4 == 14 ~.875,
      faminc_new_W4 == 15 ~.9375,
      faminc_new_W4 == 16 ~1,
      TRUE ~ NA),
    Ideology = case_when(# 1: liberal, 5: conservative
      ideo5_W4 == 1 ~ 0,
      ideo5_W4 == 2 ~ .25,
      ideo5_W4 == 3 ~ .5,
      ideo5_W4 == 4 ~ .75,
      ideo5_W4 == 5 ~ 1,
      TRUE ~ NA_real_),
    Party_ID = case_when(# 1: Democrat, 7: Republican
      pid7_W4 == 1 ~ 0,
      pid7_W4 == 2 ~ .1666667,
      pid7_W4 == 3 ~ .3333333,
      pid7_W4 == 4 ~ .5,
      pid7_W4 == 5 ~ .6666667,
      pid7_W4 == 6 ~ .8333333,
      pid7_W4 == 7 ~ 1,
      TRUE ~ NA_real_),
    Age = (2025-birthyr_W4)/100,
    Trump_Therm = Q62_2_W4/100,
    Republican_Therm = Q34_6_W4/100,
    FIRE1 = case_when(Q53_4_W4==1~0,
                      Q53_4_W4==2~.25,
                      Q53_4_W4==3~.5,
                      Q53_4_W4==4~.75,
                      Q53_4_W4==5~1),
    FIRE2 = case_when(Q53_5_W4==1~1,
                      Q53_5_W4==2~.75,
                      Q53_5_W4==3~.5,
                      Q53_5_W4==4~.25,
                      Q53_5_W4==5~0),
    FIRE3 = case_when(Q53_6_W4==1~0,
                      Q53_6_W4==2~.25,
                      Q53_6_W4==3~.5,
                      Q53_6_W4==4~.75,
                      Q53_6_W4==5~1),
    FIRE_Scale = (FIRE1+FIRE2+FIRE3)/3)

wave4 <- wave4 %>%
  mutate(ABC = if_else(is.na(ABC), 0, ABC),
         NBC = if_else(is.na(NBC), 0, NBC),
         CBS = if_else(is.na(CBS), 0, CBS),
         Fox_News = if_else(is.na(Fox_News), 0, Fox_News),
         CNN = if_else(is.na(CNN), 0, CNN),
         MSNBC = if_else(is.na(MSNBC), 0, MSNBC), 
         PBS = if_else(is.na(PBS), 0, PBS),
         Local_News = if_else(is.na(Local_News), 0, Local_News),
         Spanish_News = if_else(is.na(Spanish_News), 0, Spanish_News),
         Newsmax = if_else(is.na(Newsmax), 0, Newsmax), 
         No_News = if_else(is.na(No_News), 0, No_News))

#add recoded variables for continuous Fox Measure
wave4 <- wave4 %>%
  mutate(Fox1 = case_when(Q22_m_W4_1 ==1 ~ 1,
                          Q22_m_W4_1 ==2 ~ 0,
                          TRUE ~ 0),
         Fox2 = case_when(Q22_m_W4_2 ==1 ~ 1,
                          Q22_m_W4_2 ==2 ~ 0,
                          TRUE ~ 0),
         Fox3 = case_when(Q22_m_W4_3 ==1 ~ 1,
                          Q22_m_W4_3 ==2 ~ 0,
                          TRUE ~ 0),
         Fox4 = case_when(Q22_m_W4_4 ==1 ~ 1,
                          Q22_m_W4_4 ==2 ~ 0,
                          TRUE ~ 0),
         Fox5 = case_when(Q22_m_W4_5 ==1 ~ 1,
                          Q22_m_W4_5 ==2 ~ 0,
                          TRUE ~ 0),
         Fox6 = case_when(Q22_m_W4_6 ==1 ~ 1,
                          Q22_m_W4_6 ==2 ~ 0,
                          TRUE ~ 0),
         Fox7 = case_when(Q22_m_W4_7 ==1 ~ 1,
                          Q22_m_W4_7 ==2 ~ 0,
                          TRUE ~ 0),
         Fox8 = case_when(Q22_m_W4_8 ==1 ~ 1,
                          Q22_m_W4_8 ==2 ~ 0,
                          TRUE ~ 0),
         Fox9 = case_when(Q22_m_W4_9 ==1 ~ 1,
                          Q22_m_W4_9 ==2 ~ 0,
                          TRUE ~ 0),
         Fox10 = case_when(Q22_m_W4_10 ==1 ~ 1,
                           Q22_m_W4_10 ==2 ~ 0,
                           TRUE ~ 0),
         Fox_Total = (Fox1+Fox2+Fox3+Fox4+Fox5+Fox6+Fox7+Fox8+Fox9+Fox10)/10)

#add data on political violence
wave4 <- wave4 %>%
  mutate(Violence1 = case_when(Q64b_1_W4==1~1,
                               Q64b_1_W4==2~0.75,
                               Q64b_1_W4==3~0.5,
                               Q64b_1_W4==4~0.25,
                               Q64b_1_W4==5~0,
                               Q64b_1_W4==6~NA_real_),
         Violence2 = case_when(Q64b_2_W4==1~1,
                               Q64b_2_W4==2~0.75,
                               Q64b_2_W4==3~0.5,
                               Q64b_2_W4==4~0.25,
                               Q64b_2_W4==5~0,
                               Q64b_2_W4==6~NA_real_),
         Violence3 = case_when(Q64b_3_W4==1~1,
                               Q64b_3_W4==2~0.75,
                               Q64b_3_W4==3~0.5,
                               Q64b_3_W4==4~0.25,
                               Q64b_3_W4==5~0,
                               Q64b_3_W4==6~NA_real_))

#add data on group attitudes
wave4 <- wave4 %>%
  mutate(Latino_Therm = Q34_1_W4/100,
         White_Therm = Q34_3_W4/100,
         White_Preference = White_Therm - Latino_Therm)

#Scale evaluation: Wave 3

wave3_whites <- wave3 %>%
  filter(White==1)

grt_alpha <- wave3_whites %>%
  dplyr::select(GRT1, GRT2, GRT3)

alpha(grt_alpha)
#among whites, alpha of 0.874

fa.parallel(grt_alpha, fa="fa")
#among whites, single factor of GRT items

FIRE_alpha <- wave3_whites %>%
  dplyr::select(FIRE1,FIRE2,FIRE3)

alpha(FIRE_alpha)
#among whites, alpha of 0.617

fa.parallel(FIRE_alpha, fa="fa")
#among whites, single factor of FIRE Scale

#Confirmatory factor analysis (Wave3)

model <- 'New_GRT =~ GRT1+GRT2+GRT3
          New_FIRE =~ FIRE1+FIRE2+FIRE3'

cfa_result <- cfa(model, data = wave3_whites)

summary(cfa_result, ci = TRUE, fit.measures = TRUE)

#Correlation of GRT, FIRE, and Outgroup Prejudice among Whites: Wave 3

cor(wave3_whites$GRT_Beliefs, wave3_whites$FIRE_Scale)
#0.5906566

cor(wave3_whites$GRT_Beliefs, wave3_whites$White_Preference, use = "complete.obs")
#0.4486763

#Scale evaluation: Wave 4

wave4_whites <- wave4 %>%
  filter(White==1)

grt_alpha_4 <- wave4_whites %>%
  dplyr::select(GRT1, GRT2, GRT3)

alpha(grt_alpha_4)
#among whites, alpha of 0.876

fa.parallel(grt_alpha_4, fa="fa")
#among whites, single factor of GRT items

FIRE_alpha_4 <- wave4_whites %>%
  dplyr::select(FIRE1,FIRE2,FIRE3)

alpha(FIRE_alpha_4)
#among whites, alpha of 0.673

fa.parallel(FIRE_alpha_4, fa="fa")
#among whites, single factor of FIRE Scale

#Confirmatory factor analysis (Wave4)

model_4 <- 'New_GRT =~ GRT1+GRT2+GRT3
          New_FIRE =~ FIRE1+FIRE2+FIRE3'

cfa_result_4 <- cfa(model_4, data = wave4_whites)

summary(cfa_result_4, ci = TRUE, fit.measures = TRUE)


#Correlation of GRT, FIRE, and Outgroup Prejudice among Whites: Wave 4

cor(wave4_whites$GRT_Beliefs, wave4_whites$FIRE_Scale)
#0.6466529

cor(wave4_whites$GRT_Beliefs, wave4_whites$White_Preference, use = "complete.obs")
#0.5218506

#cross-sectional models of Fox exposure and GRT support

#dichotomous measure, wave 3
GRT_Fox_W3 <- lm(GRT_Beliefs ~ Fox_News + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites, weights = weight_white_w3)

#continuous measure, wave 3
GRT_Fox_W3_C <- lm(GRT_Beliefs ~ Fox_Total + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites, weights = weight_white_w3)

#dichotomous measure, wave 4
GRT_Fox_W4 <- lm(GRT_Beliefs ~ Fox_News + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites, weights = weight_white_p_W4)

#continuous measure, wave 4
GRT_Fox_W4_C <- lm(GRT_Beliefs ~ Fox_Total + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites, weights = weight_white_p_W4)

wordreg(l=list(GRT_Fox_W3, GRT_Fox_W3_C, GRT_Fox_W4, GRT_Fox_W4_C),
        file = "grt-models.docx",
        stars = c(0.001, 0.01, 0.05), bold = 0.05,
        custom.model.names = c("Model 1", "Model 2", "Model 3", "Model 4"),
        caption = "Predictors of Support for the GRT",
        digits = 2)

#cross-sectional analysis of waves 3-4 (no weights)

#dichotomous measure, wave 3
GRT_Fox_W3_nw <- lm(GRT_Beliefs ~ Fox_News + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites)

#continuous measure, wave 3
GRT_Fox_W3_C_nw <- lm(GRT_Beliefs ~ Fox_Total + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites)

#dichotomous measure, wave 4
GRT_Fox_W4_nw <- lm(GRT_Beliefs ~ Fox_News + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites)

#continuous measure, wave 4
GRT_Fox_W4_C_nw <- lm(GRT_Beliefs ~ Fox_Total + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites)

wordreg(l=list(GRT_Fox_W3_nw, GRT_Fox_W3_C_nw, GRT_Fox_W4_nw, GRT_Fox_W4_C_nw),
        file = "grt-models_nw.docx",
        stars = c(0.001, 0.01, 0.05), bold = 0.05,
        custom.model.names = c("Model 1", "Model 2", "Model 3", "Model 4"),
        caption = "Predictors of Support for the GRT",
        digits = 2)

#Standard Fox News Model with Each GRT Item Individually (wave 3)

GRT_Fox_grt1 <- lm(GRT1 ~ Fox_News + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites, weights = weight_white_w3)

GRT_Fox_grt2 <- lm(GRT2 ~ Fox_News + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites, weights = weight_white_w3)

GRT_Fox_grt3 <- lm(GRT3 ~ Fox_News + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites, weights = weight_white_w3)

wordreg(l=list(GRT_Fox_grt1, GRT_Fox_grt2, GRT_Fox_grt3),
        file = "grt-models_individual_items.docx",
        stars = c(0.001, 0.01, 0.05), bold = 0.05,
        custom.model.names = c("Immigrants invade and colonize the United States.", "Native-born Americans are losing their economic, political, and cultural influence in this country because of the growing population of immigrants.", "There are people who secretly work to make sure immigrants will eventually replace real Americans."),
        caption = "OLS Regression Models of Support for Individual GRT Items among White Adults, 2024 AMPS",
        digits = 2)

GRT_Fox_grt1_c <- lm(GRT1 ~ Fox_Total + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites, weights = weight_white_w3)

GRT_Fox_grt2_c <- lm(GRT2 ~ Fox_Total + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites, weights = weight_white_w3)

GRT_Fox_grt3_c <- lm(GRT3 ~ Fox_Total + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites, weights = weight_white_w3)

wordreg(l=list(GRT_Fox_grt1_c, GRT_Fox_grt2_c, GRT_Fox_grt3_c),
        file = "grt-models_individual_items_continuous.docx",
        stars = c(0.001, 0.01, 0.05), bold = 0.05,
        custom.model.names = c("Immigrants invade and colonize the United States.", "Native-born Americans are losing their economic, political, and cultural influence in this country because of the growing population of immigrants.", "There are people who secretly work to make sure immigrants will eventually replace real Americans."),
        caption = "OLS Regression Models of Support for Individual GRT Items among White Adults, 2024 AMPS",
        digits = 2)


#Standard Fox News Model with Each GRT Item Individually (wave 4)

GRT_Fox_grt1_w4 <- lm(GRT1 ~ Fox_News + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites, weights = weight_white_p_W4)

GRT_Fox_grt2_w4 <- lm(GRT2 ~ Fox_News + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites, weights = weight_white_p_W4)

GRT_Fox_grt3_w4 <- lm(GRT3 ~ Fox_News + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites, weights = weight_white_p_W4)

wordreg(l=list(GRT_Fox_grt1_w4, GRT_Fox_grt2_w4, GRT_Fox_grt3_w4),
        file = "grt-models_individual_items_w4.docx",
        stars = c(0.001, 0.01, 0.05), bold = 0.05,
        custom.model.names = c("Immigrants invade and colonize the United States.", "Native-born Americans are losing their economic, political, and cultural influence in this country because of the growing population of immigrants.", "There are people who secretly work to make sure immigrants will eventually replace real Americans."),
        caption = "OLS Regression Models of Support for Individual GRT Items among White Adults, 2024 AMPS",
        digits = 2)

GRT_Fox_grt1_w4_c <- lm(GRT1 ~ Fox_Total + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites, weights = weight_white_p_W4)

GRT_Fox_grt2_w4_c <- lm(GRT2 ~ Fox_Total + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites, weights = weight_white_p_W4)

GRT_Fox_grt3_w4_c <- lm(GRT3 ~ Fox_Total + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites, weights = weight_white_p_W4)

wordreg(l=list(GRT_Fox_grt1_w4_c, GRT_Fox_grt2_w4_c, GRT_Fox_grt3_w4_c),
        file = "grt-models_individual_items_w4_continuous.docx",
        stars = c(0.001, 0.01, 0.05), bold = 0.05,
        custom.model.names = c("Immigrants invade and colonize the United States.", "Native-born Americans are losing their economic, political, and cultural influence in this country because of the growing population of immigrants.", "There are people who secretly work to make sure immigrants will eventually replace real Americans."),
        caption = "OLS Regression Models of Support for Individual GRT Items among White Adults, 2024 AMPS",
        digits = 2)

#Hierarchically regularized entropy balancing model (simple Fox News model), wave 3

X = c("Male", "Age", "Education", "Income", "Ideology", "Party_ID", "FIRE_Scale")

hbal_w3 <- hbal(wave3_whites, "Fox_News", X, Y = "GRT_Beliefs", X.expand = X, expand.degree = 3, cv = TRUE, folds = 4, seed=12345)

summary(att(hbal_w3, method="lm_robust", dr=TRUE, displayAll=TRUE))

#note that this estimate is estimate is the sample average treatment effect on the treated (SATT), NOT the sample average treatment effect (SATE). 

#Hierarchically regularized entropy balancing model (simple Fox News model), wave 4
X = c("Male", "Age", "Education", "Income", "Ideology", "Party_ID", "FIRE_Scale")

hbal_w4 <- hbal(wave4_whites, "Fox_News", X, Y = "GRT_Beliefs", X.expand = X, expand.degree = 3, cv = TRUE, folds = 4, seed=12345)

summary(att(hbal_w4, method="lm_robust", dr=TRUE, displayAll=TRUE))

#Sensitivity Analysis, Waves 3 and 4

#Dichotomous Fox variable, Wave 3
GRT_Fox_W3 <- lm(GRT_Beliefs ~ Fox_News + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites, weights = weight_white_w3)

#Dichotomous Fox variable, Wave 4
GRT_Fox_W4 <- lm(GRT_Beliefs ~ Fox_News + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites, weights = weight_white_p_W4)

#Model to find predictors of Fox News consumption, Wave 3
Fox_Model_W3 <- glm(Fox_News ~ Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites, weights = weight_white_w3, family = "gaussian")

summary(Fox_Model_W3)

#Sensitivity analysis, Wave 3
GRT_Fox_W3_sensitivity <- sensemakr(model = GRT_Fox_W3, treatment = "Fox_News", benchmark_covariates = "FIRE_Scale", kd = 1:2, ky = 1:2, q = 1, alpha = 0.05, reduce = TRUE)

summary(GRT_Fox_W3_sensitivity)

#Model to find predictors of Fox News consumption, Wave 4
Fox_Model_W4 <- glm(Fox_News ~ Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites, weights = weight_white_p_W4, family = "gaussian")

summary(Fox_Model_W4)

#Sensitivity analysis, Wave 4
GRT_Fox_W4_sensitivity <- sensemakr(model = GRT_Fox_W4, treatment = "Fox_News", benchmark_covariates = "FIRE_Scale", kd = 1:2, ky = 1:2, q = 1, alpha = 0.05, reduce = TRUE)

summary(GRT_Fox_W4_sensitivity)

#Sensitivity Results, Wave 3

#Sensitivity Statistics:
#  Partial R2 of treatment with outcome: 0.0305 
#Robustness Value, q = 1: 0.1622 
#Robustness Value, q = 1, alpha = 0.05: 0.1034 

#Verbal interpretation of sensitivity statistics:
  
#  -- Partial R2 of the treatment with the outcome: an extreme confounder (orthogonal to the covariates) that explains 100% of the residual variance of the outcome, would need to explain at least 3.05% of the residual variance of the treatment to fully account for the observed estimated effect.

#-- Robustness Value, q = 1: unobserved confounders (orthogonal to the covariates) that explain more than 16.22% of the residual variance of both the treatment and the outcome are strong enough to bring the point estimate to 0 (a bias of 100% of the original estimate). Conversely, unobserved confounders that do not explain more than 16.22% of the residual variance of both the treatment and the outcome are not strong enough to bring the point estimate to 0.

#-- Robustness Value, q = 1, alpha = 0.05: unobserved confounders (orthogonal to the covariates) that explain more than 10.34% of the residual variance of both the treatment and the outcome are strong enough to bring the estimate to a range where it is no longer 'statistically different' from 0 (a bias of 100% of the original estimate), at the significance level of alpha = 0.05. Conversely, unobserved confounders that do not explain more than 10.34% of the residual variance of both the treatment and the outcome are not strong enough to bring the estimate to a range where it is no longer 'statistically different' from 0, at the significance level of alpha = 0.05.

#- Bounds on the Strength of Confounding: An unobserved confounder 2 time(s) as strong as the covariate FIRE_Scale could explain at most 0.99% of the residual variance of the treatment and 25.22% of the residual variance of the outcome.This means that the maximum bias that could be caused by such a confounder would bring the estimated effect to 0.0722 and the t-value to 4.216.

#Note: the Bounds on the Strength of Confounding were calculated using the coefficient estimate and se. estimate for Fox_News in the GRT_Fox_W3 model; while the t-values for the FIRE_Scale were calculated from the GRT_Fox_W3 model and Fox_Model_W3 model. The Bounds on the Strength of Confounding were estimated using the Sensemakr Shiny App, available at https://carloscinelli.shinyapps.io/robustness_value/

#Sensitivity Results, Wave 4

#Sensitivity Statistics:
#  Partial R2 of treatment with outcome: 0.0173 
#Robustness Value, q = 1: 0.1241 
#Robustness Value, q = 1, alpha = 0.05: 0.0642 

#Verbal interpretation of sensitivity statistics:
  
#  -- Partial R2 of the treatment with the outcome: an extreme confounder (orthogonal to the covariates) that explains 100% of the residual variance of the outcome, would need to explain at least 1.73% of the residual variance of the treatment to fully account for the observed estimated effect.

#-- Robustness Value, q = 1: unobserved confounders (orthogonal to the covariates) that explain more than 12.41% of the residual variance of both the treatment and the outcome are strong enough to bring the point estimate to 0 (a bias of 100% of the original estimate). Conversely, unobserved confounders that do not explain more than 12.41% of the residual variance of both the treatment and the outcome are not strong enough to bring the point estimate to 0.

#-- Robustness Value, q = 1, alpha = 0.05: unobserved confounders (orthogonal to the covariates) that explain more than 6.42% of the residual variance of both the treatment and the outcome are strong enough to bring the estimate to a range where it is no longer 'statistically different' from 0 (a bias of 100% of the original estimate), at the significance level of alpha = 0.05. Conversely, unobserved confounders that do not explain more than 6.42% of the residual variance of both the treatment and the outcome are not strong enough to bring the estimate to a range where it is no longer 'statistically different' from 0, at the significance level of alpha = 0.05.

#- Bounds on the Strength of Confounding: An unobserved confounder 2 time(s) as strong as the covariate FIRE_Scale could explain at most 0.69% of the residual variance of the treatment and 32.65% of the residual variance of the outcome.This means that the maximum bias that could be caused by such a confounder would bring the estimated effect to 0.0476 and the t-value to 3.0549.

#Note: the Bounds on the Strength of Confounding were calculated using the coefficient estimate and se. estimate for Fox_News in the GRT_Fox_W4 model; while the t-values for the FIRE_Scale were calculated from the GRT_Fox_W4 model and Fox_Model_W4 model. The Bounds on the Strength of Confounding were estimated using the Sensemakr Shiny App, available at https://carloscinelli.shinyapps.io/robustness_value/

#Write files to csv for MRobust analysis in Stata

write.csv(wave3_whites, "wave3_whites.csv")

write.csv(wave4_whites, "wave4_whites.csv")

#Importance of GRT measures
#GRT measures predict support for political violence (waves 3 and 4)
Violence1_Model_w3 <- lm(Violence1 ~ GRT_Beliefs + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites, weights = weight_white_w3)

Violence2_Model_w3 <- lm(Violence2 ~ GRT_Beliefs + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites, weights = weight_white_w3)

Violence3_Model_w3 <- lm(Violence3 ~ GRT_Beliefs + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites, weights = weight_white_w3)

Violence1_Model_w4 <- lm(Violence1 ~ GRT_Beliefs + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites, weights = weight_white_p_W4)

Violence2_Model_w4 <- lm(Violence2 ~ GRT_Beliefs + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites, weights = weight_white_p_W4)

Violence3_Model_w4 <- lm(Violence3 ~ GRT_Beliefs + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites, weights = weight_white_p_W4)

wordreg(l=list(Violence1_Model_w3, Violence2_Model_w3, Violence3_Model_w3, Violence1_Model_w4, Violence2_Model_w4, Violence3_Model_w4),
        file = "grt-models_violence.docx",
        stars = c(0.001, 0.01, 0.05), bold = 0.05,
        custom.model.names = c("Physical Harm Justified (Wave 3)", "Social Change by Any Means Necessary (Wave 3)", "Personally Willing to Cause Physical Harm (Wave 3)", "Physical Harm Justified (Wave 4)", "Social Change by Any Means Necessary (Wave 4)", "Personally Willing to Cause Physical Harm (Wave 4)"),
        caption = "Predictors of Support for Violence",
        digits = 2)

#GRT measures predict support for Trump and GOP (waves 3 and 4)

TrumpThermModel_w3 <- lm(Trump_Therm ~ GRT_Beliefs + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave3_whites, weights = weight_white_w3)

TrumpThermModel_w4 <- lm(Trump_Therm ~ GRT_Beliefs + Male + Age + Education + Income + Ideology + Party_ID + FIRE_Scale, data=wave4_whites, weights = weight_white_p_W4)

wordreg(l=list(TrumpThermModel_w3, TrumpThermModel_w4),
        file = "grt-models_Trump.docx",
        stars = c(0.001, 0.01, 0.05), bold = 0.05,
        custom.model.names = c("Warmth Toward Trump (Wave 3)", "Warmth Toward Trump (Wave 4)"),
        digits = 2)

